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ABSTRACT 

An  algorithm  is  presented  for  computing  In  z  with  complex  arithmetic, 
by  extending  to  the  complex  plane  Carlson's  treatment  of  a  classical  iteration 
using  arithmetic  and  geometric  means.  Although  not  competitive  with  current 
techniques  which  handle  the  real  and  imaginary  parts  separately,  the  algorithm 
may  be  useful  in  special  purpose  applications.  A  detailed  analysis  of 
convergence,  scaling,  and  roundoff  is  given.  Standard  identities  and  some 
minor  bookkeeping  allow  the  evaluation  of  inverse  circular  and  inverse 
hyperbolic  functions.  It  is  also  shown  that  the  basic  procedure  is  related 
to  certain  real  algorithms. 
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SIGNIFICANCE  AND  EXPLANATION 


Complex  loqarithms  are  usually  computed  with  real  arithmetic, 

12  2  -1 
ln(x  +  iy)  =  -  ln(x  +  y  )  +  i  tan  (y/x)  , 

using  the  real  In  and  real  tan  *  functions  available  in  high  level  programming 
languages.  This  report  studies  a  novel  approach,  based  on  complex  arithmetic, 
obtained  by  extending  to  the  complex  plane  a  known  real  algorithm.  The  strategy 
consists  of  a  basic  iteration,  which  uses  one  complex  square  root  per  cycle, 
accelerated  by  Richardson  extrapolation,  a  speeding-up  process  well-known  to 
numerical  analysts.  With  proper  scaling  of  the  arqument,  the  accelerated 
procedure  requires  four  complex  square  roots  for  accuracy  to  10  decimals 
(compared  to  eighteen  complex  square  roots  in  the  non-acoelerated  case) . 

Standard  identities  such  as 

-1  f~2 

cosh  z=ln(z+/z  -1)  , 

-1  .  "1 
cos  z  =  i  cosh  z  . 

etc . 

and  some  minor  bookkeeping  to  account  for  principal  branches,  allow  the  evalua¬ 
tion  of  complex  inverse  trigonometric  functions.  The  report  analyzes  the 
algorithm  with  respect  to  scaling,  convergence,  and  stability,  and  it  relates 
the  algorithm  to  other  procedures,  including  classical  methods  for  calculating 
the  constant  r, .  Since  complex  square  roots  are  done  at  the  software  level,  and 
thus  costly  in  machiri'*  time,  the  algorithm  is  not  competitive  with  standard 
methods  handling  the  real  and  imaginary  parts  separately  with  real  arithmetic. 
However,  its  simplicity  and  stability  could  make  it  attractive  for  implementation 
m  microcode  or  road  only  memory  in  special  purfiose  applications  requiring 
extensive  u  .in  of  •  •  1  nm,  inter  y  cornel  ex  functions.  - - -  ... 


EVALUATION  OF  COMPLEX  LOGARITHMS  AND  RELATED  FUNCTIONS 

T 

George  J.  Miel 

1.  Introduction.  Logarithms  of  complex  numbers  are  commonly  computed 
using  real  arithmetic  separately  for  the  real  and  imaginary  parts, 

12  2  -1 

ln(x  +  iy)  =  —  In (x  +  y  )  +  i  tan  (y/x)  ,  (1.11 

with  suitable  precautions  to  avoid  numerical  problems;  see,  e.g.,  [C,  Algorithm 
243].  We  analyze  an  algorithm  based  on  complex  arithmetic,  obtained  by 
extending  to  the  complex  case  Carlson's  procedure  [3].  The  strategy  consists 
of  a  basic  iteration,  which  uses  one  complex  square  root  per  cycle,  accelerated 
by  Richardson  extrapolation.  The  basic  iteration  generates  one  of  the  sequences 
of  Borchardt’s  algorithm  [2,  p.  499],  [4,  p.  170).  For  real  arguments,  this 
iteration  is  also  related  to  Thacher's  algorithm  for  inverse  cosines  [11],  to 
Viete's  infinite  product  for  it  [8,  p.  26],  and  to  the  method  of  equal  perim¬ 
eters  [8,  p.  32] .  As  for  the  real  case,  the  improvement  due  to  extrapolation 
is  substantial,  the  algorithm  is  reliable  and  stable,  and  storage  needs  are 
modest  as  there  are  no  constants  to  be  saved.  Numerical  experiments  indicate 
no  serious  cancellation  leading  to  loss  of  significant  figures,  as  sometimes 
happens  when  a  real  algorithm  is  extended  to  the  complex  case.  With  an  adequate 
reduction  in  the  range  of  the  independent  variable,  the  accelerated  procedure 
requires  four  complex  square  roots  for  lOD  accuracy.  Standard  identities,  and 
some  bookkeeping  to  account  for  principal  branches,  allow  the  evaluation  ::f 
inverse  circular  and  inverse  hyperbolic  functions.  The  complex  arithmetic 
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provides  a  unified  design  for  a  simple  and  modular  software  package.  Unfortu¬ 
nately,  the  complex  square  roots  preclude  the  algorithm  from  competition  with 
the  straightforward  approach  (1.1),  which  takes  advantage  of  efficient  real 
elementary  functions  provided  in  high  level  languages.  However,  the  algorithm 
might  be  of  interest  in  special  purpose  applications  implemented  in  microcode 
or  read  only  memory.  Certain  monotonicity  properties  of  Borchardt's  iteration 
can  be  exploited  for  computation  in  complex  interval  arithmetic  [10] . 
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2.  The  Basic  Algorithm.  The  goal  is  to  evaluate  the  single-valued  function 
In  z  whose  range  is  {x  +  iy  |  -  tt  <  y  <  it}.  The  rational  operations  of  complex 
arithmetic  and  the  principal  complex  square  root  are  assumed  available  lb,  V t.y 
The  latter  is  the  single-valued  function  whose  range  is 

C+  =  (x  +  iy  |  x  >  0  or  x=0  with  y  S  0}  . 


Recall  that 


f~T  .  r  + 

/z  =  z,  if  z  e  C 


(2.1) 


We  will  use  the  function  F(z)  =  (z  coth  z)/w  and  its  expansion 

.  .  2  4  2m  „ ,  2m  +  2, 

F(z)  =  aQ  +  a  z  +  a^z  +  ...  +  a^z  +  0(z  ), 


(2.2) 


where  a  ^  depends  on  the  Bernoulli  number  ^ , 


22^B  . 

a  =  - 2J_  , 

3  (2j>:  w 

and  where  Landau's  notation  f(z)  =  0(g(z))  is  adequately  defined  [7,  p.  156). 
In  what  follows,  D  =  {z  |  |z|  >  l} . 


Theorem.  Let  z  e  D  and 


z+1 


i-l  '  ’n+1  sn 


=  g  +  /f2-  1,  n  ^  0 . 


(2.1) 


The  numbers  b  satisfy 
n 


„-n- 1 

=  coth  2  w. 


=  2^  +  f.  , 

n+1  n  n 


(2.4) 


where  w  =  In  z,  =  -1/g^  ,  I  r'n  I 


b  I ,  and  the  numbers  u  =  2  n  r 
n  n  '>n 


converge  to  u  =  1/w  with  a  rate 


u-u  ,  =-  4_1  (u-u  )  +  0(16  n  ’). 
n  +  1  n 


(2.5) 
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Proof.  For  n  =  0,  5  =  (eW  +  l)/(eW-  1)  =  coth  w/2.  The  relation 

o 

__  2 

=  a  +  ib,  csch  £  =  |  sinh  £  |  (sinh  a  cos  b  -  i  cosh  a  sin  b)  implies  that 

w ,'2  a  C+  whenever  w  e  ln(D)  .  Use  the  identities  coth  £/2  =  coth  £  + 

,  2  2 

cscn  ,  csch  £  =  coth  £  -  1  and  (2.1)  to  get 

coth  w/4  =  coth  w/2  +  / coth2  w/2  -  1,  w  £  ln(D) . 

The  relation  (2.3)  gives  5  ,  =  25  -  1/5  The  convergence  follows  from 

n+1  n  n+1 

u  =  F ( 2  n  'w)  and  lim  F(£)  =  1/w.  Finally,  use  (2.2)  with  m  =  1  to  get  (2.5) 
n  £**0 

The  second  expression  in  (2.4)  shows  that  as  n  increases,  5  ,  gets 

n+1 

increasingly  close  to  25  .  This  fact  provides  a  simple  variable  precision 

n 

scheme.  A  range  reduction  allows  the  evaluation  of  logarithms  for  machine 
representable  arguments. 

Basic  Algorithm.  Adequate  precision  is  assumed  available.  Given  z'  /  0, 
proceed  as  follows  to  find  In  z'  correct  to  d  decimals: 

1.  Factorize  |z'|  =  2rx  ,  x  e  [2,4).  Let  z  =  2  ' . 

2.  Compute  (2.3)  with  n  =  0,  1,  ...,  N  where  N  is  such  that  5..  and 

N 

25..  .  agree  to  d  +  1  decimals. 

N-i 

3.  Let  In  z'  =  (2"N_15„)“1  +  r  In  2. 

N 

Various  range  reductions  are  possible.  Instead  of  the  modulus,  one  can 
use  | a  +  ib |  =  |a|  +  |b|  or  | a  +  ib |  =  max  (|a|,  |b|).  The  code  should  take 
advantage  of  the  multiplications  by  powers  of  2. 

Roundoff  Propagation.  Numerical  experiments  indicate  that  the  algorithm 
is  remarkably  stable.  A  simplified  analysis  shows  why.  Consider 


5*5  -  e  , 

o  o  o 


'  .  =  <M[  )  -  5  , 

n+1  ^  n  n 


-n-lc 

u  =  2  5  , 

n  n 


where  cf)  ( z )  =  z  +  /z2  -  1  and  6  reflects  the  accuracy  of  the  complex  square 

n 

A 

root  routine.  Letting  e  =  £  -  £  > 

n  n  n 

£  4>  ’  (C  )E  +  <s  »  <t>'(£  )  =  5  ,/(£  1  -  £  ). 

n+1  n  n  n  n  n+1  n+1  n 

Since  £  ,  =  2£  ,  assume  for  simplicity  that  6  =26  ,  ,  d> '  (£  )  =  2 . 

n+1  n  n  n-1  n 

Then  e  =  N2N-16  +  2N£  ,  and  assuming  no  error  in  the  multiplication  by 

N  o  o 

-N-l 

2  ,  we  get 

~  n  -  en 

U  "U  *  TO  +  °  . 

N  N  4  o  2 

The  accumulated  roundoff  is  acceptable  if  the  square  root  routine  is  accurate 
to  at  least  d  +  2  decimals.  Practice  shows  that  cancellation  causes  actual 


roundoff  to  be  smaller  than  (2.6). 
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3.  Richardson  Extrapolation.  Section  5  shows  that  {u^}  is  one  of  the 

sequences  generated  by  Borchardt's  algorithm  C 2] ,  [3],  [4,  p.  170] .  As  for 

the  real  case,  (2.5)  indicates  that  successive  errors  in  u  are  ultimately 

n 

reduced  by  a  factor  1/4.  In  order  to  speed-up  convergence,  we  extend  to  the 

complex  case  the  treatment  of  Carlson  C  3] . 

The  procedure  is  given  by 

P  =u,  0  :£  n  £  N, 

on  n 


P  =  P  +  - 

kn  k-l,n+l 


’k-l.n+l  ~  pk-l,n  1  <  k  <  N 

4k  _  x  ’  0  <  n  <  N-k  * 


The  scheme  generates  a  triangular  array , 


P00  y  p01  P02  "•  ?0,N-1  /  P0n 

iXX  IX 

T)  T3  ^  T7 


P  P 

.10  ^  11 


1,N-1 


P  P 

,20  2 ,  N-  2 


in  which  the  arrows  illustrate  the  dependence  of  P  on  P  and  P 

Kn  k-i  ,  n  k- i , n+i 


Recall  that  for  z  e  D,  u  =  F  (£  )  u  *  1/w, 

n  n 


=  w/2n+1. 


Fq (C )  =  F (C) ,  Fk(;) 


=  4*Fk-l<0  -  Fk.1(2g 


4-1 


Then  (2.2)  implies  that 


Fk(g  =  U  +  bkX+2  +  0(;2k+4) , 


where  b,  is  a  constant.  Since  ?  =  2  £  ,  ,  we  have  P,  =  F ,  (C  ,  )  .  Consequently, 

k  n  n+1  kn  k  n+k 

P  -  „  ,  4-<k+1)  lk+”*1,b  «2k+2  *  0(4-(k+2>  (k™+1>  ). 

kn  k 

Thus,  the  errors  in  successive  elements  of  the  k-th  row  of  (3.3)  are  each  time 
roughly  reduced  by  a  factor  4  „ 

Error  Bounds  via  Interpolation  Theory.  Let  f(z)  =  F  (/"z)  =  (/"z  coth  /z ) /w 
and  z^  =  £2  =  w2/4n+\-  (3.2)  is  the  Neville  scheme  for  evaluating  polynomials 

with 

P  =  value  at  z=0  of  the  polynomial  of  degree  £  k 
kn 


which  interpolates  f(z)  at  z  ,  z  z  ,  ; 

n  n+1  n+k 

see  Brezinski  Cl,  p.26].  The  function  f(z)  is  analytic  for  |z|  < 


Theorem.  If  u  =  1/w,  w  =  In  z,  and 

/l-n1  -  18  'h 

e  f 


1  <  z  S  e 


(3.4) 


then 


|u  _  P  I  <  2“(k+1)  (k+2n+1)  |w|2k+1/M' 

1  kn 1 


+  tt]  . 


(3.5) 


Proof.  Let  a  =  and  let  C"  denote  the  square  with  vertices 
a(±l±i).  We  first  show  that 


max  | coth  z\  ^  coth  a. 
zeC" 


(3.6) 


For  the  side  z  =  a  +  iy,  -a  ^  y  5  a,  we  get 


coth  z  < 


lea+ly[  +  |e~a~lyl 

[|ea+iy|  _  |e-a-iY| 

ia  1 , 


a  -a 
e  +  e 


a  -a 

e  -  e 


=  coth  a. 


For  z  =  x  +  ia,  -a  <  x  <  a,  use  e  =  — (-1  +  i  / 3)  to  obtain 
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Estimate  max  |  /t  coth  / 1 J  for  t  r  C'  rather  thari  t  t  C  usina  that 
/t  is  then  on  the  square  C". 


| co th  z 

Similarly  for  the  other  two  sides.  The  function  'Jj(x)  reaches  its  maximum  at 
x=±a.  Since  iH±a)  <  coth  a,  we  get  (3.6). 

Now,  let  C  denote  the  union  of  two  semicircles  of  radius  1/2  and  centers 

2 

w 

at  0  and  zQ  =  —  respectively,  and  two  parallel  segments  as  shown  in  the 
figure.  Consider  also  the  union  of  two  partial  parabolas, 

C'  =  {±  (s2-  a2  +  2ias)|  -a  £  s  Z  a}. 

The  curve  C'  is  inside  the  circle  of  analyticity  and  (3.4)  implies  that  C 
is  not  outside  C' . 

By  a  classic  result. 


1  -  i/3  tanh  x 

( 3  tanh  x  +  1 

tanh  x  -  i  /3 

t  tanh  x  +  3  ■ 

u 


'kn 


(-1) 


k+1 


z  z 
n  n+1 


n+k 


2tt  i 


w 


I, 


( 


where 


I  =/ 


g(t)  dt 


ct(t-zn)...(t-zn+k) 


g(t)  =  /t  coth  /IT. 


But 


ill  £  2k+2(2|zQ|  +  IT)  M, 
where  M  k  pax  [ g(t)|  .  We  have 

£fg  |g(t)  I  <  ^a*,  |g(t)  I  =  ma^„  |g(z2)  j  <  /  2  a  coth  a  «=  0.9  72  f, 

where  we  used  (3.6).  To  complete  the  proof,  take  M  =  tt  and  use  (3.7)  and  v  3 . ... 
For  the  range  reduction  described  earlier,  (3.4)  is  satisfied  and 
w !  <  |  In  4  +  i  tt  |.  For  accuracy  to  d  decimals,  set  the  bound  in  (3.5)  no 
greater  than  —10  ^  and  thus  obtain: 
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Non-accelerated  case  (k=0,n=N)  N^  1.66  d  +  2.48, 

Accelerated  case  (k=n  =  N)  N^  /l.ll  d  +  1.66  -  0.07. 

Actual  computation  shows  that  10D  accuracy  requires  N  =  18  in  the  non-accel¬ 
erated  case  and  N  =  4  in  the  extrapolated  case.  The  acceleration  defined 
by  (3.1)  and  (3.2)  is  easily  coded.  The  code  should  take  advantage  of  the 
multiplications  by  powers  of  4.  Variable  precision  is  possible  by 
a  priori  specification  of  N  or  by  on-line  comparison  of  P  and 


P 

N-1,1 

Roundoff  Propagation.  If  the  values  in  the  first  row  of  (3.3)  are 
contaminated  with  errors  whose  magnitudes  are  less  than  e ,  then  the  errors 
later  in  the  extrapolation  have  magnitude  which  nowhere  exceed  2e .  Combining 
this  with  (2.6),  we  get 


P  -p 

1  NO  NO ! 


<  ii<50  +  g  ♦  id 


where  P  „  is  the  computed  value  of  P„,0  and  C  is  the  accumulated  roundoff 
NO  Nu 

in  the  computation  of  (3.2).  The  major  term  is  and  the  extrapolation 

is  well-conditioned  provided  that  the  square  root  routine  is  of  good  quality. 


1 
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4,  Inverse  Trigonometric  Functions.  The  table  below  specifies  for  each 
function  a  bijection  between  the  domain  and  given  range. 


FUNCTION 

RANGE 

In 

R1  = 

{a+ib  |  -IT  S>  b  < 

t} 

cosh  1 

R2  = 

{a+ib  |  -tt  <  b  < 

0  or  b  =  -t,  0  with  a  £  ol 

sinh 

II 

m 
C X 

{a+ib  |  -7T/2  <  b 

<  t/2  or  b  =  ±  t/2  with  a  > 

tanh  "*■ 

R4  = 

R1  -  {—  it  /  2 } 

-1 

cos 

1R2 

.  -1 

sin 

iR3 

-1 

tan 

-iR4 

The  following  procedures  use  standard  identities  to  evaluate  the  functions. 


Algorithm  for  w  =  cosh  1z. 


Algorithm  for  w  =  sinh  z. 


w'=ln(z+/z  -1)  =a+ib. 
If  b>0  then  w=-w ' . 


w '  =ln  (z+/z2  +1)  =a+ib . 

If  be  [  -t,-t/2)  then  w=-iiT-w' 
If  b=-t,  0  and  a<0  then  w=-a+ib.  If  be  (7T/ 2 , tt )  then  w=it-w'. 


Otherwise,  w=w' . 


Algorithm  for  w  =  cos  ~*~z 


If  b=±t/2  and  a<0  then  w=-a+ib. 
Otherwise,  w=w'. 

Algorithm  for  w  =  sin  "''z. 


. -1 

w=  i  cosh  z. 


w=  i  sinh  (-iz) 


-1 


Algorithm  for  w  =  tanh  z. 

z?±\,  w=”  In  (j— —■)  • 

2  1-  Z 


-1 


Algorithm  for  w  =  tan  z. 
z/±i,  w=-itanh  1iz. 


In  the  case  of  cosh  ^  and  sinh  \  some  logic  is  needed  in  order  to  choose 
the  value  in  the  specified  range,  since  these  functions  are  double-valued 
in  the  range  of  In,  F  9,  p.  417]. 
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5,  Connection  with  Real  Algorithms.  For  real  arguments,  our  accelerated 

v 

procedure  is  exactly  equivalent  to  Carlson's  treatment  of  Borchardt's  algorithm 

The  non-accelerated  procedure  is  related  to  Thacher's  algorithm  [11]  for  real 
inverse  cosines  and  to  classical  methods  for  calculating  tt . 


Borchardt ' s  Algorithm .  If  u  ,  v  ^  >0  then  the  sequences 


u  ..  =  -r(u  +v  )  , 
n+1  2  n  n 


v  =  /u  v  ,  n  2  -1, 
n+1  n+1  n 


(5.1) 


converge  monotonically  to  a  common  limit  B(u_^,v  ^)  ,  [2].  We  show  that 

Lr  _  P  /TS  2  2~ 

5n'  Cn+1  =  ?n+  /Sn  -  U-l+  V-1  * 


u  =  2 
n 


-n-1. 


(5.2) 


Get  the  invariant  4n+1(u2  ,  -  v2  , )  =  c  and  then  substitute  v  =  (u2  -  4  nc)1//2 

n+1  n+1  n  n 

in  the  first  relation  of  (5.1).  Our  basic  algorithm  generates  (5.2)  corres¬ 


ponding  to  B 


z2+l 


2z 


2  _  2  . 
z  -1  z  -1 


In  z 


Thacher's  Algorithm.  If  R  =  / 2z+2,  R  =  / R  +2  then  t  =  2IV [ R  —  2 | 
-  — 2 -  n+l  n  n  n 

converge  to  cos  ^z  if  |z|  <  1  and  to  cosh  ^z  if  |z|  SI,  [11].  It  turns 

out  that  t^  is  the  reciprocal  of  v^  generated  by  (5.1)  with  proper  u_^  and  v 


Method  of  Equal  Perimeters 

B 


If  V 


1  ^  TT 

—  cot  — , 

U  U 


>  2  then 


1  it 

—  esc  — 

U  y 


u  =  ■  cot  — =  radius  of  inscribed  circle  in  a  regular 

n  2n+iy  2n+  y 

2n+^p  -  gon  of  perimeter  2. 


Descartes  worked  with  U  =  4,  [8,  p.32]. 


1 


Viete' s  Infinite  Product.  If  g_^ 
then  lim(g_1gQ. .  .g  )  =  4/tt  ,  [8,  p.26]. 
corresponds  to  B(l/  /  2)  =  1/sin  ^"]f  ' 


We  have 


-,/l  1  7T 

l  9  +  1  5 


g-iV 


•  g„  * 


v  ,  wher 
n 
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